﻿// 提取点数据

#include <windows.h>
#include <iostream>
#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include <stdexcept>
#include <string>
#include <sstream>
#include <fstream>
#define UNICODE
#define _UNICODE
#include "tinydir.h"
#include <cstdlib>
#include "wdatetime.h"
using namespace std;
using std::cout;

#include "blh2ecef.h"
using namespace blh2ecef;

#include <netcdf.h> // for atmPrf
#include "wstringutils.h"





//二维散点值
struct ScatterPointValue {
    double _x = 0;
    double _y = 0;
    double _alt = 0;
    double _value = 0;
};

bool extractPointsByNC(string filename, vector<ScatterPointValue>& resultScatters, string& error)
{
    int ncid = 0 , varid = 0;

    int data_in[111] ;

    int retval = nc_open(filename.c_str(), NC_NOWRITE, &ncid) ;
    if (retval != 0) {
        cout << "Error code from nc_open " << retval << endl;
        return false;
    }
    
    int lonId = 0;
    int latId = 0;
    int altId = 0;
    int tempId = 0;
    retval = nc_inq_varid(ncid, "Lon", &lonId);
    if (retval != 0) {
        cout << "failed to get variable Lon" << endl;
        return false;
    }
    retval = nc_inq_varid(ncid, "Lat", &latId);
    if (retval != 0) {
        cout << "failed to get variable Lat" << endl;
        return false;
    }
    retval = nc_inq_varid(ncid, "MSL_alt", &altId );
    if (retval != 0) {
        cout << "failed to get variable MSL_alt" << endl;
        return false;
    }
    retval = nc_inq_varid(ncid, "Temp", &tempId);
    if (retval != 0) {
        cout << "failed to get variable Temp" << endl;
        return false;
    }

    int ndimsp = 0;
    int nvarsp = 0;
    int nattsp = 0;
    int unlimdimidp = 0;

    retval = nc_inq(ncid,
        &ndimsp,
        &nvarsp,
        &nattsp,
        &unlimdimidp);

    int ndims[1] = { 0 };
    int dimids[1] = { 0 };
    //nc_inq_dimids(ncid, ndims, dimids, 0);
    size_t shape[1] = { 0 };
    nc_inq_dimlen(ncid, dimids[0],shape);
    if (shape[0] == 0) {
        cout << "bad shape value of dimension 0" << endl;
        return false;
    }
    vector<float> londata(shape[0]);
    vector<float> latdata(shape[0]);
    vector<float> altdata(shape[0]);
    vector<float> tempdata(shape[0]);

    nc_get_var_float(ncid, lonId, londata.data() );
    nc_get_var_float(ncid, latId, latdata.data());
    nc_get_var_float(ncid, altId, altdata.data());
    nc_get_var_float(ncid, tempId, tempdata.data());
    nc_close(ncid);

    for (int i = 0; i < tempdata.size(); ++i) {
        ScatterPointValue s;
        s._alt = altdata[i]*1000.0;
        s._x = londata[i];
        s._y = latdata[i];
        s._value = tempdata[i];
        resultScatters.push_back(s);
    }

    return true;
}

bool extractPoints(bool useIon, string filename,vector<ScatterPointValue>& resultScatters , string& error)
{//对一个剖面文件
    GDALDataset* ipds = (GDALDataset*)GDALOpenEx(filename.c_str(), GDAL_OF_VECTOR, NULL, NULL, NULL);
    if (ipds == 0) {
        error = "bad prf file.";
        return false;
    }
    OGRLayer* lyr = ipds->GetLayer(0);
    string lyrName = lyr->GetName();
    int index1 = 0;
  
    for (auto& feature : lyr) {
        OGRGeometry* poGeometry = feature->GetGeometryRef();
        OGRPoint* pt = poGeometry->toPoint();
        //int nf = feature->GetFieldCount();
        //for (int iif = 0; iif < nf; ++iif) {
        //    OGRFieldDefn* defn = feature->GetFieldDefnRef(iif);
        //    string name = defn->GetNameRef();
        //    cout << iif<<" field name " << name << endl;
        //}
        double val = 0.0;
        if (useIon)val = feature->GetFieldAsDouble("TEC_cal"); //TEC_cal  TEC_rel  ELEC_dens
        else val = feature->GetFieldAsDouble("Temp"); 
        //printf("%d %.3f %.3f %.3f %.3f \n", index1, pt->getX(), pt->getY(), pt->getZ(), ion);//lon(deg) , lat(deg) , altitude(km)
        index1++;
        {
            double z = pt->getZ()*1000.0;
            ScatterPointValue spv;
            spv._x = pt->getX();
            spv._y = pt->getY();
            spv._alt = z;
            spv._value = val;
            resultScatters.push_back(spv);
        }
    }
    GDALClose(ipds);
    return true;
}

double xyzDist2(blh2ecef::CoordX left, blh2ecef::CoordX right) {
    double x = left.x - right.x;
    double y = left.y - right.y;
    double z = left.z - right.z;
    double dist = x * x + y * y + z * z;
    return dist;
}

int main(int argc , char* argv[] )
{
    std::cout << "version v0.3" << endl;
    cout << "A program to extract points from ion file. 2025/7/8" << endl;
    cout << "usage:yy_ionprofPoints i|a mindistkm /some/input/dir filenamePrefix outname.geojson" << endl;
    cout << " mindistkm is min distance(km) between two points, smaller distance will get rid of." << endl;
    cout << " a is for atmPrf, i is for ionPrf " << endl;
    cout << "提取通过文件名前缀控制获取哪些文件的廓线数据。" << endl;
    cout << "ionPrf的高程数据单位是KM。" << endl;
    GDALAllRegister();
    if (argc != 6) {
        cout << "Error , argc not 6. argc=" << argc << endl;
        return 11;
    }

    string prfType = argv[1];//  "a";
    double mindistkm = atof(argv[2]) * 1000 ;
    const double MIN_DIST2 = mindistkm * mindistkm;//两个点最小距离，少于这个距离跳过。
    // atm  E:/data/YunYao/ZhanTing/atmprf20250704/20250704
    // ion  E:/data/YunYao/ZhanTing/ionPrf/ionPrf
    string indir = argv[3];// "E:/data/YunYao/ZhanTing/atmprf20250704/20250704";
    string inprefix = argv[4];//  "atmPrf_Y005";// "ionPrf_Y005";
    string outname = argv[5];//  "y005_atm_output.geojson";
    string ecef_outname = wStringUtils::replaceString( outname , ".geojson", "_ecef.geojson") ;

    
    bool useIon = true;
    if (prfType.compare("a") == 0) useIon = false;//use Atm

    //{"dt":"2025-07-15T12:00:00Z",
    //"lon":100.29826,
    //"lat":40.96040,
    //"alt":  1022.000
    //
    //cesium
    //"x":-862441.8699218825,
    //"y":4746530.447596349,
    //"z":4159773.1618734878 检查过了 C++ blh2ecef 计算结果和cesium一样得。
    ////blh2ecef::CoordX testXYZ = blh2ecef::blh2ecef({ 40.96040, 100.29826 , 1022.000 });

    vector<string> filenames;
    vector<string> satNames;
    vector<string> baseNames;
    tinydir_dir dir;
    wchar_t wdir[512] = {};
    size_t newSize = 0;

    mbstowcs_s(&newSize, wdir, indir.size() + 1, indir.c_str(), _TRUNCATE);
    tinydir_open(&dir, wdir);

    while (dir.has_next)
    {
        tinydir_file file;
        tinydir_readfile(&dir, &file);

        //wchar to char
        size_t rlen = 0;
        char cName[512] = {};
        wcstombs_s(&rlen, cName, (size_t)511, file.name, (size_t)511);
        // printf("%s\n", cName);// cName is baseName no dir path.
        string bName = string(cName);
        string fullName = indir + "/" + bName;

        if (file.is_dir == false && bName.find(inprefix) != string::npos)
        {
            filenames.push_back(fullName);
            string satN = bName.substr(7, 4);
            satNames.push_back(satN);
            baseNames.push_back(bName);
        }
        tinydir_next(&dir);
    }
    tinydir_close(&dir);


    string jsonText = "{\"type\":\"FeatureCollection\",\"name\":\"temp-multi-lines\",\"crs\":"
        "{\"type\":\"name\",\"properties\":{\"name\":\"urn:ogc:def:crs:OGC:1.3:CRS84\"}},\"features\":[";

    string jsonText_ecef = "{\"type\":\"FeatureCollection\",\"name\":\"temp-multi-lines\",\"crs\":"
        "{\"type\":\"name\",\"properties\":{\"name\":\"ECEF\"}},\"features\":[";

    int lineCnt = 0;
    int totalPtCount = 0;
    for (int ifile = 0; ifile < filenames.size(); ++ifile) {
        //一条线
        string bname = baseNames[ifile];
        int year = atof(bname.substr(12, 4).c_str());
        int dayOfYear = atof(bname.substr(17, 3).c_str());
        int hour = atof(bname.substr(21, 3).c_str());
        int minu = atof(bname.substr(24, 3).c_str());

        int ymd8 = wDatetime::convertYmd7toYmd8(year*1000+dayOfYear);
        int month = (ymd8 % 10000) / 100;
        int day = ymd8 % 100;

        char dtstr[64];
        sprintf(dtstr, "%04d-%02d-%02dT%02d:%02d:00Z", year, month, day, hour, minu);
        string dtstr2(dtstr);

        vector<ScatterPointValue> scatters;
        string error;
        bool ok = false;
        if (useIon) {
            ok = extractPoints(useIon, filenames[ifile],
                scatters, error);
        }
        else {
            ok = extractPointsByNC(filenames[ifile], scatters, error);
        }
        if (scatters.size() > 1) { //一个点以上的线才有效
            stringstream ss, ss_ecef;
            ss << "{\"type\":\"Feature\",\"properties\":{\"sat\":\""
                << satNames[ifile]
                <<"\""
                <<",\"dt\":\""<< dtstr2 <<"\"" 
                <<"},\"geometry\":{\"type\":\"LineString\",\"coordinates\":[";
            ss_ecef << ss.str();
            blh2ecef::CoordX ecef0 = { 0,0,0 };
            int newPtCount = 0;
            for (int is = 0; is < scatters.size(); ++is) {
                //点太多这里要精简到1km一个点吧
                blh2ecef::CoordX ecefv = blh2ecef::blh2ecef({ scatters[is]._y, scatters[is]._x , scatters[is]._alt });
                double dist2 = xyzDist2(ecef0, ecefv);
                if (dist2 > MIN_DIST2) {
                    if (is > 0) {
                        ss << ",";
                        ss_ecef << ",";
                    }
                    ss << "[" << scatters[is]._x << "," << scatters[is]._y << "," << scatters[is]._alt << "," << scatters[is]._value << "]";
                    ss_ecef << "[" << ecefv.x << "," << ecefv.y << "," << ecefv.z << "," << scatters[is]._value << "]";
                    ecef0 = ecefv;
                    newPtCount++;
                }
            }
            ss << "]}}";
            ss_ecef << "]}}";
            if (newPtCount > 1) {
                //超过一个点才有效
                if (lineCnt > 0) {
                    jsonText += ",";
                    jsonText_ecef += ",";
                }
                jsonText += ss.str();
                jsonText_ecef += ss_ecef.str();
                lineCnt++;
                if (lineCnt % 10 == 0) {
                    cout << ", lineCnt " << lineCnt << std::flush;
                }
                totalPtCount += newPtCount;
            }

        }
    }

    jsonText+="]}";
    jsonText_ecef += "]}";
    ofstream ofs(outname.c_str());
    ofs << jsonText;
    ofs.close();

    ofstream ofs_ecef(ecef_outname.c_str());
    ofs_ecef << jsonText_ecef;
    ofs_ecef.close();

    cout << "total vertex count " << totalPtCount <<" ,  line count "<<lineCnt  << endl;
    if (lineCnt > 0) {
        cout << "average pts per line : " << totalPtCount / lineCnt << endl;
    }

    cout << "done" << endl;

    return 0;
}
